/*=======================================================================================
	Heterogeneity.do
	
		Uses the Personal Finance Website (PFW) data to analyze the response of 
		consumption to the APFD payments across different subpopulations
		
	Author: Lorenz Kueng, June 2015
=========================================================================================*/

cap log close PFW_07a
log using "$homedir/log-files/PFW_07a_$date.log", text replace name(PFW_07a)


*=============================================================
* Programs
*=============================================================

cap n program drop MakeTable
program MakeTable

	cap n test APFD_q1 = APFD_q5
		local Pval_15 = `r(p)'
		local Fval_15 = `r(F)'
	cap n test APFD_q1 = APFD_q4
		local Pval_14 = `r(p)'
		local Fval_14 = `r(F)'
	cap n test APFD_q1 = APFD_q3
		local Pval_13 = `r(p)'
		local Fval_13 = `r(F)'
	cap n test APFD_q1 = APFD_q2
		local Pval_12 = `r(p)'
		local Fval_12 = `r(F)'
	cap n test APFD_q2 = APFD_q5				
		local Pval_25 = `r(p)'
		local Fval_25 = `r(F)'
	cap n test APFD_q2 = APFD_q4				
		local Pval_24 = `r(p)'
		local Fval_24 = `r(F)'
	cap n test APFD_q2 = APFD_q3
		local Pval_23 = `r(p)'
		local Fval_23 = `r(F)'
	cap n test APFD_q3 = APFD_q5
		local Pval_35 = `r(p)'
		local Fval_35 = `r(F)'
	cap n test APFD_q3 = APFD_q4
		local Pval_34 = `r(p)'
		local Fval_34 = `r(F)'
	cap n test APFD_q4 = APFD_q5
		local Pval_45 = `r(p)'
		local Fval_45 = `r(F)'

	forvalues i=1(1)5 {		
		cap lab var APFD_q`i' "PFD payments x quantile `i'"
	}

	outreg2 using "${file}.xls",  alpha(0.01, 0.05, 0.1) symbol(***,**,*) bdec(3) se nocons label  addstat(Number of clusters,e(N_clust)) ///
		ctitle("${name}") ///
		addtext( ///
				" - Time FE", Yes, /// 
				" - Family size FE", Yes, ///
				" - State FE", Yes, /// 
				" - Quantile FE", Yes, ///
				" - Income change", Yes, ///
				" - Household characteristics", Yes, ///
				" - F-value 1=5", "`Fval_15'", ///						
				" - P-value 1=5", "`Pval_15'", ///						
				" - F-value 1=4", "`Fval_14'", ///						
				" - P-value 1=4", "`Pval_14'", ///						
				" - F-value 1=3", "`Fval_13'", ///						
				" - P-value 1=3", "`Pval_13'", ///						
				" - F-value 1=2", "`Fval_12'", ///						
				" - P-value 1=2", "`Pval_12'", ///						
				" - F-value 2=5", "`Fval_25'", ///						
				" - P-value 2=5", "`Pval_25'", ///						
				" - F-value 2=4", "`Fval_24'", ///						
				" - P-value 2=4", "`Pval_24'", ///						
				" - F-value 2=3", "`Fval_23'", ///						
				" - P-value 2=3", "`Pval_23'", ///						
				" - F-value 3=5", "`Fval_35'", ///						
				" - P-value 3=5", "`Pval_35'", ///						
				" - F-value 3=4", "`Fval_34'", ///						
				" - P-value 3=4", "`Pval_34'", ///						
				" - F-value 4=5", "`Fval_45'", ///						
				" - P-value 4=5", "`Pval_45'" ///						
		) ///
		keep( APFD_q* )
end



*=============================================================
* Data
*=============================================================

use "$homedir/data/stata/PFW_quarterly.dta", clear

	generate date = Quarter
	format   date %tq

	lab var APFD "PFD payments" // APFD payments indentified via direct deposit entry in early October

	xtset userid date


*** Create additional variables

	** Controls
	foreach var in residentialstatus marriagestatus {
		display _n(1) "`var'"
		cap drop `var'_
		encode   `var', generate(`var'_)
		bysort userid (date): egen _temp = mode(`var'_), maxmode
		replace  `var'_ = _temp if `var'_==. & _temp!=. 
		replace  `var'_ = 999   if `var'_==. 
		drop _temp
	}
	foreach var in age edu profession homezipcode {
		display _n(1) "`var'"
		cap drop `var'_
		generate `var'_ = `var'
		bysort userid (date): egen _temp = mode(`var'_), maxmode
		replace  `var'_ = _temp if `var'_==. & _temp!=. 
		replace  `var'_ = 999   if `var'_==. 
		drop _temp
	}

	cap drop _temp
	gen D_incomeQ = D.incomeQ
	bysort userid (Quarter): egen _temp = mean(D_incomeQ)
	replace D_incomeQ = _temp if missing(D_incomeQ)==1
	drop _temp*

	foreach var in incomeQ D_incomeQ incomeY liquid_asset1 {
		winsor `var', g(`var'W) p(0.01)
		sum `var' `var'W, d
	}

	** Total expenditures as a measure of permanent income (excluding uncategorized transactions, ie only nondurables, services, durables)

		bysort userid: egen totexp_annual_nominal = mean(totexp1_nominal) // average quarterly total spending (nominal)
		bysort userid: egen totexp_annual = mean(totexp1) // average quarterly total spending (real)
		replace totexp_annual_nominal = -totexp_annual_nominal*4 // annualize and make expenditures positive
		replace totexp_annual = -totexp_annual*4
		winsor  totexp_annual_nominal, p(0.01) generate(totexp_annual_nominalW)
		winsor  totexp_annual, p(0.01) generate(totexp_annualW)



*===============================================================================
* MPC Heterogeneity - Table  4: "MPC Heterogeneity"
*					- Figure 6: "Welfare Loss"
*===============================================================================

sort userid date

lab var APFD "PFD payments"

*
foreach depvar in nondur_serv {

	global depvar = "`depvar'"

	cap log close D${depvar}_heterog
	log using "$homedir/results/tables/Table4_D${depvar}_Heterogeneity_${date}.log", text replace name(D${depvar}_heterog)


		global file "$homedir/results/tables/Table4_D${depvar}_Heterogeneity_${date}"

		cap rm "${file}.xls"
		cap rm "${file}.txt"

		xtset userid date


		*-----------------------------------------------------------
		* Winsorize dependent variable at 1%
		*-----------------------------------------------------------

		* save a copy of the original variable

		qui cap drop ${depvar}_orig
		qui generate ${depvar}_orig = ${depvar}

		* winsorize at 1% (conditional on non-zero observations for disaggregated categories)

		count if ${depvar}==0 

		if `r(N)'<5000 {
		
			cap drop temp
			winsor ${depvar}, generate(temp) p(0.01) 
			replace temp = 0 if temp==.
			replace ${depvar} = temp
		}
		else { // winsorize conditional on non-zeros for disaggregated variables

			cap drop temp
			winsor ${depvar} if ${depvar}!=0, generate(temp) p(0.01) 
			replace temp = 0 if temp==.
			replace ${depvar} = temp
		}

		* make all variables positive (i.e., expenditures)

		qui replace $depvar = abs(${depvar}) 
		summarize ${depvar}, detail


		*-----------------------------------------------------------------------
		* Average effect (baseline): just to check - not for table
		*-----------------------------------------------------------------------

		*** Impose same sample in all regressions

			reghdfe ${depvar} APFD ///
					family_sizeImp ///
					liquid_asset1 ///
					incomeY ///
					D_incomeQ ///
				if APFDid_annual!=. /// NOTE: APFDid_annual!=. conditions on Alaskans for whom we know they received a dividend via direct deposit in early October
				, absorb(date userid marriagestatus_ residentialstatus_ age_ edu_ profession_ homezipcode_) vce(cluster userid)

			cap drop _e_sample
			generate _e_sample = e(sample)


		*-----------------------------------------------------------------------
		* Liquidity
		*-----------------------------------------------------------------------

		*** Liquidity (bank account balance)

			egen Liquidity_quantile = xtile(liquid_asset1), n(5) by(year Alaska) 

			** Quantiles			
			cap drop quantile
			generate quantile = Liquidity_quantile
			global name "Liquid Assets"

			** Interact with quantiles 			
			cap drop I_quantile*
			xi i.quantile, prefix(I_) noomit

			forvalues i=1(1)5 {		
				cap drop APFD_q`i'
				generate APFD_q`i' = APFD * I_quantile_`i'
			}

			** Regression 
			reghdfe D.${depvar} APFD_q* ///
					D_incomeQ /// 
				if _e_sample==1 /// 
				, absorb(date Alaska family_sizeImp age_ edu_ marriagestatus_ residentialstatus_ profession_ homezipcode_ quantile) vce(cluster userid)			

			MakeTable



		*** Relative liquidity (cash-on-hand): bank balances / quarterly total spending

			generate LiquidityRelative = liquid_asset1_nominal/(totexp_annual_nominalW/4)

			egen LiquidityRelative_quantile = xtile(LiquidityRelative), n(5) by(year Alaska) 

			** Quantiles			
			cap drop quantile
			generate quantile = LiquidityRelative_quantile
			global name "Relative Liquidty"

			** Interact with quantiles 			
			cap drop I_quantile*
			xi i.quantile, prefix(I_) noomit

			forvalues i=1(1)5 {		
				cap drop APFD_q`i'
				generate APFD_q`i' = APFD * I_quantile_`i'
			}

			** Regression 
			reghdfe D.${depvar} APFD_q* ///
					D_incomeQ /// 
				if _e_sample==1 /// 
				, absorb(date Alaska family_sizeImp age_ edu_ marriagestatus_ residentialstatus_ profession_ homezipcode_ quantile) vce(cluster userid)			

			MakeTable



		*-----------------------------------------------------------------------
		* By income
		*-----------------------------------------------------------------------

		*** Current income per capita

			** Income per capita  (household-equivalent)
			generate incomeYWpercap = incomeYW / equivalence
			bysort userid: egen _temp = mean(incomeYWpercap) // annualize
			replace incomeYWpercap = _temp
			drop _temp
			xtset userid date

			egen incomeYW_quantile = xtile(incomeYWpercap), n(5) by(year Alaska)

			** Quantiles			
			cap drop quantile
			generate quantile = incomeYW_quantile
			global name "Income"

			** Interact with quantiles 			
			cap drop I_quantile*
			xi i.quantile, prefix(I_) noomit

			forvalues i=1(1)5 {		
				cap drop APFD_q`i'
				generate APFD_q`i' = APFD * I_quantile_`i'
			}	

			** Regression 
			xtset userid date
			reghdfe D.${depvar} APFD_q* ///
					D_incomeQ /// 
				if _e_sample==1 /// 
				, absorb(date Alaska family_sizeImp age_ edu_ marriagestatus_ residentialstatus_ profession_ homezipcode_ quantile) vce(cluster userid)			

			MakeTable

			bysort quantile: summarize incomeYWpercap if Alaska==1
			
			
				** Store results for Figure
				
				xtset userid date

				cap drop _quantile
				generate _quantile =.
				 
				foreach var in  _beta _se _2se_up _2se_down {
					cap drop `var'_${name}
					generate `var'_${name} =.
				}
				
				forvalues i=1/5 {
					replace _quantile = `i'                            		   in `i'
					replace _beta_${name}     = _b[APFD_q`i']                  in `i'
					replace _se_${name}       = _se[APFD_q`i']                 in `i'
					replace _2se_up_${name}   = _b[APFD_q`i']+2*_se[APFD_q`i'] in `i'
					replace _2se_down_${name} = _b[APFD_q`i']-2*_se[APFD_q`i'] in `i'
				}
			
				lab var _beta_${name} "by current income"

				cap drop mlab*		
				generate mlabMPC_${name} = string(_beta_${name},"%3.2f")
			
		

		*** Permanent income per capita

			generate totexp_annual_percap = totexp_annual / equivalence
			egen totexpY_quantile = xtile(totexp_annual_percap), n(5) by(year Alaska) 

			** Quantiles			
			cap drop quantile
			generate quantile = totexpY_quantile
			global name "Totexp"

			** Interact with quantiles 			
			cap drop I_quantile*
			xi i.quantile, prefix(I_) noomit

			forvalues i=1(1)5 {		
				cap drop APFD_q`i'
				generate APFD_q`i' = APFD * I_quantile_`i'
			}	
			
			*xtset userid date

			** Regression 
			xtset userid date
			reghdfe D.${depvar} APFD_q* ///
					D_incomeQ /// 
				if _e_sample==1 /// 
				, absorb(date Alaska family_sizeImp age_ edu_ marriagestatus_ residentialstatus_ profession_ homezipcode_ quantile) vce(cluster userid)			

			MakeTable

			bysort quantile: summarize totexp_annual_percap  if Alaska==1
			bysort quantile: summarize incomeYWpercap        if Alaska==1


				** Store results for Figure
				
				xtset userid date
				
				foreach var in  _beta _se _2se_up _2se_down {
					cap drop `var'_${name}
					generate `var'_${name} =.
				}
				
				forvalues i=1/5 {
					replace _beta_${name}     = _b[APFD_q`i']                  in `i'
					replace _se_${name}       = _se[APFD_q`i']                 in `i'
					replace _2se_up_${name}   = _b[APFD_q`i']+2*_se[APFD_q`i'] in `i'
					replace _2se_down_${name} = _b[APFD_q`i']-2*_se[APFD_q`i'] in `i'
				}

				
				generate mlabMPC_${name} = string(_beta_${name},"%3.2f")


			*=======================================================================
			* Figure 7a: MPC by income (current and permanent)
			*=======================================================================
				
				lab var _beta_Totexp "MPC by permanent income"
				lab var _beta_Income "MPC by current income"
				
				twoway ///
					   connected _beta_Totexp          _quantile, mc(dkgreen) ms(O) lc(dkgreen) lw(medthick) connect(line) mlabel(mlabMPC_Totexp) mlabc(dkgreen) mlabposition(12) || ///
					   line      _2se_up_Totexp        _quantile, lc(dkgreen) lp(dash) lw(thin) || ///
					   line      _2se_down_Totexp      _quantile, lc(dkgreen) lp(dash) lw(thin) || ///
					   connected _beta_Income          _quantile, mc(black) ms(X)  lc(black) lw(medthick) connect(line) mlabel(mlabMPC_Income) mlabc(black) mlabposition(12) || ///
					   line      _2se_up_Income        _quantile, lc(black) lp(dash) lw(thin) || ///
					   line      _2se_down_Income      _quantile, lc(black) lp(dash) lw(thin) || ///
					, xtitle("income per capita quintiles") xlabel(1(1)5) ///
					  ylabel(0.0(0.2)0.8,angle(0) labcolor(black) tlcolor(black) axis(1)) yscale(lcolor(black) axis(1)) ytitle("MPC",color(black)) ///
					  yline(0, lc(black) lp(solid) lw(thin) axis(1)) ///
					  graphregion(color(white)) legend(order(1 4) col(1) ring(0) pos(12) ) 
					  
					local file "$homedir/results/figures/PFW_D${depvar}_Heterogeneity_by_income"
						graph export "`file'.eps", replace
						cap rm       "`file'.pdf"
						!epstopdf    "`file'.eps"
						rm           "`file'.eps"
						graph export "`file'.tif", replace


			*=======================================================================
			* Table R-2: PFD payments by income quintile for Referee 3
			*=======================================================================

				cap log close TableR2
				log using "$homedir/results/tables/TableR2_PFDpayments_by_IncomeQuintiles_${date}.log", text replace name(TableR2)

					forvalues q=1/5 {
						tabstat APFD_q`q' if APFD_q`q'>0, stat(N mean sd)
					}

				log close TableR2 
				
				
			
		*-----------------------------------------------------------------------
		* By shock size
		*-----------------------------------------------------------------------	

		*** Relative shock size (using current annualized income)

			generate ShockSizeY = APFD_annual_nominal/incomeYW
			replace  ShockSizeY = ShockSizeY * 0.76 if state=="WA" // attenuation factor 2010-2014

			egen ShockSizeY_quantile = xtile(ShockSizeY), n(5) by(year Alaska) // Note: by year and state since we have to impute PFD for users in Washington based on self-reported family size instead of number of observed PFD checks

			** Quantiles
			cap drop quantile
			generate quantile = ShockSizeY_quantile
			global name "Relative Size,by income"

			** Interact with quantiles 			
			cap drop I_quantile*
			xi i.quantile, prefix(I_) noomit
			
			forvalues i=1(1)5 {		
				cap drop APFD_q`i'
				generate APFD_q`i' = APFD * I_quantile_`i'
			}

			** Regression 
			reghdfe D.${depvar} APFD_q* ///
					D_incomeQ /// 
				if _e_sample==1 /// 
				, absorb(date Alaska family_sizeImp age_ edu_ marriagestatus_ residentialstatus_ profession_ homezipcode_ quantile) vce(cluster userid)			

			MakeTable



		*** Relative shock size (using permanent income)

			generate ShockSize = APFD_annual_nominal/totexp_annual_nominalW
			replace  ShockSize = ShockSize * 0.76 if state=="WA" // attenuation factor 2010-2014

			egen ShockSize_quantile = xtile(ShockSize), n(5) by(year Alaska) // Note: by year and state since we have to impute PFD for users in Washington based on self-reported family size instead of number of observed PFD checks

			** Quantiles
			cap drop quantile
			generate quantile = ShockSize_quantile
			global name "Relative Size,per totexp"

			** Interact with quantiles 			
			cap drop I_quantile*
			xi i.quantile, prefix(I_) noomit

			forvalues i=1(1)5 {		
				cap drop APFD_q`i'
				generate APFD_q`i' = APFD * I_quantile_`i'
			}

			** Regression 
			reghdfe D.${depvar} APFD_q* ///
					D_incomeQ /// 
				if _e_sample==1 /// 
				, absorb(date Alaska family_sizeImp age_ edu_ marriagestatus_ residentialstatus_ profession_ homezipcode_ quantile) vce(cluster userid)			

		MakeTable



		*** Dividend payment size 

			local l1 = 1000
			local l2 = 1800
			local l3 = 2500
			local l4 = 3700
			cap drop _temp
			generate _temp = .
			replace  _temp = 1 if APFD_annual_nominal>=0    & APFD_annual_nominal<`l1'
			replace  _temp = 2 if APFD_annual_nominal>=`l1' & APFD_annual_nominal<`l2'
			replace  _temp = 3 if APFD_annual_nominal>=`l2' & APFD_annual_nominal<`l3'
			replace  _temp = 4 if APFD_annual_nominal>=`l3' & APFD_annual_nominal<`l4'
			replace  _temp = 5 if APFD_annual_nominal>=`l4' & APFD_annual_nominal<.
	
			tab _temp if _e_sample==1 & Alaska==0			

			** Quantiles			
			cap drop quantile
			generate quantile = _temp
			global name "Dividend Size"

			** Interact with quantiles 			
			cap drop I_quantile*
			xi i.quantile, prefix(I_) noomit

			forvalues i=1(1)5 {		
				cap drop APFD_q`i'
				generate APFD_q`i' = APFD * I_quantile_`i'
			}

			** Regression 
			reghdfe D.${depvar} APFD_q* ///
					D_incomeQ /// 
				if _e_sample==1 /// 
				, absorb(date Alaska family_sizeImp age_ edu_ marriagestatus_ residentialstatus_ profession_ homezipcode_ quantile) vce(cluster userid)			

			MakeTable



		*** Family size 
					
			cap drop _temp
			generate _temp = 1 if family_size==1
			replace  _temp = 2 if family_size==2
			replace  _temp = 3 if family_size==3
			replace  _temp = 4 if family_size==4
			replace  _temp = 5 if family_size >4

			tab _temp if _e_sample==1			

			** Quantiles			
			cap drop quantile
			generate quantile = _temp
			global name "Family Size"

			** Interact with quantiles 			
			cap drop I_quantile*
			xi i.quantile, prefix(I_) noomit
			
			forvalues i=1(1)5 {		
				cap drop APFD_q`i'
				generate APFD_q`i' = APFD * I_quantile_`i'
			}	
			
			** Regression 
			reghdfe D.${depvar} APFD_q* ///
					D_incomeQ /// 
				if _e_sample==1 /// 
				, absorb(date Alaska age_ edu_ marriagestatus_ residentialstatus_ profession_ homezipcode_ quantile) vce(cluster userid)			
	
			MakeTable	



		*-----------------------------------------------------------------------
		* Using squared instead of relative dividend as in Scholnick (2013)
		*-----------------------------------------------------------------------

			generate APFD_squared = (APFD/100)^2
			lab var  APFD_squared "(PFD payments/100)^2"
			global name "Dividend Squared"

			reghdfe D.${depvar} APFD APFD_squared ///
					D_incomeQ /// 
				if _e_sample==1 /// 
				, absorb(date Alaska family_sizeImp age_ edu_ marriagestatus_ residentialstatus_ profession_ homezipcode_) vce(cluster userid)			

					outreg2 using "${file}.xls",  alpha(0.01, 0.05, 0.1) symbol(***,**,*) bdec(3) se nocons label  addstat(Number of clusters,e(N_clust)) ///
						ctitle("${name}") ///
						addtext( ///
								" - Time FE", Yes, /// 
								" - Family size FE", Yes, ///
								" - State FE", Yes, /// 
								" - Quantile FE", --, ///
								" - Income change", Yes, ///
								" - Household characteristics", Yes, ///
								" - Time FE", Yes, /// 
								" - Family size FE", Yes, ///
								" - State FE", Yes, /// 
								" - Income", Yes, ///
								" - Liquid assets", Yes, ///
								" - Household characteristics", Yes ///						
						) ///
					keep( APFD APFD_squared )

}

save "_temp_PFW_07a", replace





*===============================================================================
* MPC Heterogeneity - Bootstrap standard errors of ex-post loss
*===============================================================================

	*----------------------------------
	* Bootstrap SEs 
	*----------------------------------

use "_temp_PFW_07a", clear

	*** Model parameters
	
	global depvar = "nondur_serv"
	global T= 4 
	global gamma = 2
	set seed 123456789


	*** Keep only necessary variables to speed up bootstrap

			** Regression 
			reghdfe D.${depvar} APFD_q* ///
					D_incomeQ /// 
				if _e_sample==1 /// 
				, absorb(date Alaska family_sizeImp age_ edu_ marriagestatus_ residentialstatus_ profession_ homezipcode_ quantile) vce(cluster userid)			

	keep ${depvar} userid date year Alaska APFD ShockSize D_incomeQ _e_sample family_sizeImp age_ edu_ marriagestatus_ residentialstatus_ profession_ homezipcode_  // ShockSize_quantile D_incomeQ APFDid_annual APFD APFD_q* totexp1_nominal APFD_annual_nominal year state liquid_asset1

	save "_temp_bootstrap", replace

	
	** Create dataset to store bootstrap results

	cap drop PotentialLossQ*
	cap drop ActualLossQ*
	cap drop MPCQ
	forvalues q=1/5 {
		generate PotentialLossQ`q'=.
		generate ActualLossQ`q'=.
		generate MPCQ`q'=.
	}
	
	keep PotentialLossQ* ActualLossQ* MPCQ*
	drop if _n>=1
	save "$homedir/data/stata/economic_loss_bootstrap_results", replace


forvalues b=0/1000 { // number of bootstrap replications

	display _n(5)"*** bootstrap sample `b' : $S_TIME  $S_DATE ****"

	use "_temp_bootstrap", clear

	xtset userid date

	if `b'>0 {
		bsample, cluster(userid) idcluster(newid) // create a bootstrap sample
	}
	if `b'==0 {
		generate newid = userid
	}

	xtset newid date
 

	*** Potential loss

		cap drop PotentialLoss
		generate PotentialLoss = (ShockSize)^2 * (${T}-1) * (${gamma}/2)


	*** Quintiles 

		** compute relative dividend and potential-loss quintiles
		cap drop ShockSize_quantile
		egen     ShockSize_quantile = xtile(ShockSize), n(5) by(year Alaska) // Note: by year and state since we have to impute PFD for users in Washington based on self-reported family size instead of number of observed PFD checks

		cap drop PotentialLoss_quantile
		egen     PotentialLoss_quantile = xtile(PotentialLoss), n(5) by(year Alaska) 
		replace  PotentialLoss_quantile =. if Alaska==0

		** interact PFD with potential-loss quintiles 	
		cap drop quantile
		generate quantile = ShockSize_quantile
		cap drop I_quantile*
		xi i.quantile, prefix(I_) noomit
		forvalues i=1(1)5 {
			cap drop APFD_q`i'
			generate APFD_q`i' = APFD * I_quantile_`i' 
			lab var  APFD_q`i' "`i'"
		}


	*** Response by potential-loss quintile

		xtset newid date

		reghdfe D.${depvar} APFD_q* ///
				D_incomeQ /// 
			if _e_sample==1 /// 
			, absorb(date Alaska family_sizeImp age_ edu_ marriagestatus_ residentialstatus_ profession_ homezipcode_ quantile) vce(cluster newid)			


	*** Store results

		forvalues q=1/5 {

			cap drop MPCQ`q'
			generate MPCQ`q'=.
			replace  MPCQ`q'= _b[APFD_q`q'] if Alaska==1 & I_quantile_`q'==1

			cap drop PotentialLossQ`q'
			generate PotentialLossQ`q'=.
			replace  PotentialLossQ`q'= PotentialLoss if Alaska==1 & I_quantile_`q'==1

			cap drop ActualLossQ`q'
			generate ActualLossQ`q'=.
			replace  ActualLossQ`q'= PotentialLossQ`q' * (MPCQ`q')^2

			replace PotentialLossQ`q' = 100* PotentialLossQ`q' // in percent 
			replace ActualLossQ`q'    = 100* ActualLossQ`q'    // in percent 
			
			cap drop MPCQ`q'_SE0
			if `b'==0	generate MPCQ`q'_SE0 = _se[APFD_q`q'] // initial SE (clustered, not bootstrapped)
			if `b'>0	generate MPCQ`q'_SE0 = .
			
		}

	keep if Alaska==1
	keep MPCQ* PotentialLossQ* ActualLossQ* quantile
	collapse (mean) MPCQ* PotentialLossQ* ActualLossQ*, by(quantile)
	collapse (firstnm) MPCQ* PotentialLossQ* ActualLossQ*

	generate bootstrap_count = `b'
	order *_SE0
	order bootstrap_count MPCQ1-MPCQ5 *_SE0 PotentialLossQ* ActualLossQ*

	append using "$homedir/data/stata/economic_loss_bootstrap_results", force

	save "$homedir/data/stata/economic_loss_bootstrap_results", replace

} // end bootstrap	

rm "_temp_bootstrap.dta"



	*==================================
	* Figure 7b 
	*==================================

use "$homedir/data/stata/economic_loss_bootstrap_results", clear

	*** Compute 95% confidence intervals

		ds MPC* PotentialLoss* ActualLoss*
		foreach var in `r(varlist)' {
			gen `var'_p5  = `var'
			gen `var'_p95 = `var'
			gen `var'_med = `var'
		}

		sort bootstrap_count
		collapse (first) MPCQ1-MPCQ5 MPCQ*_SE0 PotentialLossQ1-PotentialLossQ5 ActualLossQ1-ActualLossQ5 (p5) *_p5 (p95) *_p95 (median) *_med
		drop MPCQ*_SE0_*

		set obs 5
					
		foreach var in MPC PotentialLoss ActualLoss {
			gen `var'     =.
			gen `var'_p5  =.
			gen `var'_p95 =.
			gen `var'_med =.

			forvalues i=1/5 {

				sum  `var'Q`i'
				replace `var' = `r(mean)' in `i'
				drop `var'Q`i'

				sum  `var'Q`i'_p5
				replace  `var'_p5 = `r(mean)' in `i'
				drop `var'Q`i'_p5

				sum  `var'Q`i'_p95
				replace  `var'_p95 = `r(mean)' in `i'
				drop `var'Q`i'_p95

				sum  `var'Q`i'_med
				replace  `var'_med = `r(mean)' in `i'
				drop `var'Q`i'_med			
			}
		}

		generate MPC_SE0 =., after(MPC)
		
		forvalues i=1/5 {

			sum  MPCQ`i'_SE0
			replace  MPC_SE0 = `r(mean)' in `i'
			drop MPCQ`i'_SE0
		}
		
		generate MPC_2SEup   = MPC + 1.96*MPC_SE0, after(MPC_SE0)
		generate MPC_2SEdown = MPC - 1.96*MPC_SE0, before(MPC_2SEup)
		
		gen quantile = _n
		order quantile, first


	*** Figures

		cap drop mlab*		
		generate mlabMPC = string(MPC,"%3.2f")
		generate mlabPot = string(PotentialLoss,"%3.2f")
		generate mlabAct = string(ActualLoss,"%3.2f")
		*replace  mlabExP = "" in 1

		lab var PotentialLoss "potential loss"
		lab var MPC "MPC by size"
		lab var ActualLoss "actual loss"

		twoway ///
			   connected MPC               quantile, mc(dkgreen) ms(O) lc(dkgreen) connect(line) mlabel(mlabMPC) mlabc(dkgreen) mlabposition(12) || ///
			   line      MPC_2SEdown       quantile, lc(dkgreen) lp(dash) || ///
			   line      MPC_2SEup         quantile, lc(dkgreen) lp(dash) || ///
			   /// line      MPC_p5            quantile, lc(dkgreen) lp(dash) || ///
			   /// line      MPC_p95           quantile, lc(dkgreen) lp(dash) || ///
			   connected PotentialLoss     quantile, mc(maroon)  ms(S) lc(maroon) connect(line) mlabel(mlabPot) mlabc(maroon) mlabposition(10) yaxis(2) || ///
			   line      PotentialLoss_p5  quantile, lc(maroon)  lp(dash) yaxis(2) || ///
			   line      PotentialLoss_p95 quantile, lc(maroon)  lp(dash) yaxis(2) || ///
			   connected ActualLoss        quantile, mc(blue)    ms(X) lc(blue) connect(line) mlabel(mlabAct) mlabc(blue) mlabposition(6) yaxis(2) || ///
			   line      ActualLoss_p5     quantile, lc(blue)    lp(dash) yaxis(2) || ///
			   line      ActualLoss_p95    quantile, lc(blue)    lp(dash) yaxis(2) || ///
			, xtitle("relative dividend size quintiles") xlabel(1(1)5) ///
			  ylabel(0.0(0.2)1.2,angle(0) labcolor(dkgreen) tlcolor(dkgreen) axis(1)) yscale(lcolor(dkgreen) axis(1)) ytitle("MPC",color(dkgreen)) ///
			  /// yline(0, lc(dkgreen) lp(line) lw(thin) axis(1)) ///
			  ylabel(0.0(1.0)5.0,angle(0) labcolor(blue)    tlcolor(blue)    axis(2)) yscale(lcolor(blue)    axis(2)) ytitle("Economic Loss (in %)",color(blue) orientation(rvertical) axis(2)) ///
			  /// yline(0, lc(blue) lp(line) lw(thin) axis(2)) ///
			  graphregion(color(white)) legend(order(1 4 - 7) col(2) ring(0) pos(12) )

			local file "$homedir/results/figures/MPC_by_PotentialLoss"
				graph export "`file'.eps", replace
				cap rm       "`file'.pdf"
				!epstopdf    "`file'.eps"
				rm           "`file'.eps"
				graph export "`file'.tif", replace


		*** Figures for slides

		lab var PotentialLoss "potential loss"
		//lab var MPC "MPC by size"
		lab var MPC "   "
		//lab var ActualLoss "actual loss"
		lab var ActualLoss "           "

		twoway ///
			   connected MPC               quantile, mc(none) ms(none)  lc(none) connect(none) /*mlabel(mlabMPC)*/ mlabc(none) mlabposition(12) || ///
			   line      MPC_2SEdown       quantile, lc(none) lp(dash) || ///
			   line      MPC_2SEup         quantile, lc(none) lp(dash) || ///
			   connected PotentialLoss     quantile, mc(maroon)  ms(S)  lc(maroon) connect(line) mlabel(mlabPot) mlabc(maroon) mlabposition(10) yaxis(2) || ///
			   line      PotentialLoss_p5  quantile, lc(maroon)  lp(dash) yaxis(2) || ///
			   line      PotentialLoss_p95 quantile, lc(maroon)  lp(dash) yaxis(2) || ///
			   connected ActualLoss        quantile, mc(none)    ms(none)  lc(none) connect(none) /*mlabel(mlabAct)*/ mlabc(none) mlabposition(6) yaxis(2) || ///
			   line      ActualLoss_p5     quantile, lc(none)    lp(dash) yaxis(2) || ///
			   line      ActualLoss_p95    quantile, lc(none)    lp(dash) yaxis(2) || ///
			, xtitle("relative dividend size quintiles") xlabel(1(1)5) ///
			  ylabel(0.0(0.2)1.2,angle(0) labcolor(dkgreen) tlcolor(dkgreen) axis(1)) yscale(lcolor(dkgreen) axis(1)) ytitle("MPC",color(dkgreen)) ///
			  ylabel(0.0(1.0)5.0,angle(0) labcolor(blue)    tlcolor(blue)    axis(2)) yscale(lcolor(blue)    axis(2)) ytitle("Economic Loss (in %)",color(blue) orientation(rvertical) axis(2)) ///
			  graphregion(color(white)) legend(order(1 4 - 7) col(2) ring(0) pos(12) )

			local file "$homedir/results/figures/forSlides/MPC_by_PotentialLoss_1"
				graph export "`file'.eps", replace
				cap rm       "`file'.pdf"
				!epstopdf    "`file'.eps"
				rm           "`file'.eps" 


		lab var PotentialLoss "potential loss"
		lab var MPC "MPC by size"
		//lab var ActualLoss "actual loss"
		lab var ActualLoss "           "

		twoway ///
			   connected MPC               quantile, mc(dkgreen) ms(O)  lc(dkgreen) connect(dkgreen) mlabel(mlabMPC) mlabc(dkgreen) mlabposition(12) || ///
			   line      MPC_2SEdown       quantile, lc(dkgreen) lp(dash) || ///
			   line      MPC_2SEup         quantile, lc(dkgreen) lp(dash) || ///
			   connected PotentialLoss     quantile, mc(maroon)  ms(S)  lc(maroon) connect(line) mlabel(mlabPot) mlabc(maroon) mlabposition(10) yaxis(2) || ///
			   line      PotentialLoss_p5  quantile, lc(maroon)  lp(dash) yaxis(2) || ///
			   line      PotentialLoss_p95 quantile, lc(maroon)  lp(dash) yaxis(2) || ///
			   connected ActualLoss        quantile, mc(none)    ms(none)  lc(none) connect(none) /*mlabel(mlabAct)*/ mlabc(none) mlabposition(6) yaxis(2) || ///
			   line      ActualLoss_p5     quantile, lc(none)    lp(dash) yaxis(2) || ///
			   line      ActualLoss_p95    quantile, lc(none)    lp(dash) yaxis(2) || ///
			, xtitle("relative dividend size quintiles") xlabel(1(1)5) ///
			  ylabel(0.0(0.2)1.2,angle(0) labcolor(dkgreen) tlcolor(dkgreen) axis(1)) yscale(lcolor(dkgreen) axis(1)) ytitle("MPC",color(dkgreen)) ///
			  ylabel(0.0(1.0)5.0,angle(0) labcolor(blue)    tlcolor(blue)    axis(2)) yscale(lcolor(blue)    axis(2)) ytitle("Economic Loss (in %)",color(blue) orientation(rvertical) axis(2)) ///
			  graphregion(color(white)) legend(order(1 4 - 7) col(2) ring(0) pos(12) )

			local file "$homedir/results/figures/forSlides/MPC_by_PotentialLoss_2"
				graph export "`file'.eps", replace
				cap rm       "`file'.pdf"
				!epstopdf    "`file'.eps"
				rm           "`file'.eps"


log close PFW07a
